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Abstract 

Uncertainty  quantification  has  been  an  active  research  area  in  the  past  15  years  be¬ 
cause  of  its  potential  of  significant  applications  ranging  from  signal  processing  to  aircraft 
wing  designs.  It  is  well  understood  that  effective  numerical  methods  for  stochastic  par¬ 
tial  differential  equations  (SPDES)  are  essential  for  uncertainty  quantification.  In  the 
last  decade  much  progress  has  been  made  in  the  construction  of  numerical  algorithms 
to  efficiently  solve  SPDES  with  random  coefficients  and  white  noise  perturbations. 
However,  high  dimensionality  and  low  regularity  are  still  the  bottleneck  in  solving  real 
world  applicable  SPDES  with  efficient  numerical  methods.  This  project  is  intended  to 
address  the  numerical  analysis  as  well  as  algorithm  aspects  of  SPDES.  Three  major 
contributions  are  made  in  this  project:  i)  Construction  and  convergence  analysis  of 
Quasi  Monte  Carlo  based  Particle  Swarm  Optimization  (PSO)  method;  ii)  Efficient 
adaptive  domain  sparse  grid  method  for  SPDES;  iii)  High  order  methods  of  SPDES 
via  systems  of  forward  backward  stochastic  differential  equations.  Our  work  contains 
algorithm  constructions,  rigorous  error  analysis,  and  extensive  numerical  experiments 
to  demonstrate  our  algorithm  efficiency  and  validity  of  our  theoretical  analysis. 


1  Introduction 

The  overall  goal  of  our  AFOSR  sponsored  research  program  is  to  construct  fast  and  efficient 
numerical  algorithms  for  solving  stochastic  partial  differential  equations  and  apply  them  to 
solve  nonlinear  filtering  problems  and  optimal  control  problems  under  SPDE  constraints. 

Uncertainty  quantification  has  been  an  active  research  area  in  the  past  15  years  with  many 
significant  application  potentials  ranging  from  signal  processing  to  aircraft  wing  designs.  It 
is  well  understood  that  effective  numerical  methods  for  stochastic  partial  differential  equa¬ 
tions  (SPDES)  are  essential  for  uncertainty  quantification.  In  the  last  decade  much  progress 
has  been  made  in  the  construction  of  numerical  algorithms  to  efficiently  solve  SPDES  with 
random  coefficients  and  white  noise  perturbations.  However,  high  dimensionality  and  low 
regularity  are  still  the  bottleneck  in  solving  real  world  applicable  SPDES  with  efficient  numer¬ 
ical  methods.  This  project  is  intended  to  address  the  numerical  analysis  as  well  as  algorithm 
aspects  of  SPDES.  The  following  are  our  main  contributions  during  the  funding  period. 

(i).  Quasi  Monte  Carlo  based  Particle  Swarm  Optimization  (PSO)  method.  Inspired 
by  the  social  behavior  of  the  bird  flocking  or  fish  schooling,  the  particle  swarm  optimization 
(PSO)  is  a  population  based  stochastic  optimization  method  developed  by  Eberhart  and 
Kennedy  in  1995.  It  has  been  used  across  a  wide  range  of  applications.  Faure,  Halton  and 
Vander  Corput  sequences  have  been  used  for  initializing  the  swarm  in  PSO.  Quasirandom  (or 
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low-discrepancy)  sequences  such  as  Faure,  Halton,  Vander  Corput  etc  are  deterministic  and 
suffers  from  correlations  between  radical  inverse  functions  with  different  bases  used  for  dif¬ 
ferent  dimensions.  In  this  paper,  we  investigate  the  effect  of  initializing  the  swarm  with 
scrambled  optimal  Halton  sequence,  which  is  a  randomized  quasirandom  sequence.  This  en¬ 
sures  that  we  still  have  the  uniformity  properties  of  quasirandom  sequences  while  preserving 
the  stochastic  behavior  for  particles  in  the  swarm.  Numerical  experiments  are  conducted 
with  benchmark  objective  functions  with  high  dimensions  to  verify  the  convergence  and 
effectiveness  of  the  proposed  initialization  of  PSO. 

(ii) .  Efficient  adaptive  domain  sparse  grid  method  for  SPDES  The  key  of  constructing 
efficient  numerical  algorithms  for  SPDES  is  to  overcome  the  obstacles  of  low  regularity  and 
high  dimensionality,  and  the  associated  unbounded  domain  problem.  In  our  construction,  we 
used  the  importance  sampling  method  to  construct  bounded  domains  for  the  Zakai  equation 
with  a  very  small  number  of  sample  points.  Then  we  used  the  split  up  scheme  to  construct 
higher  order  numerical  algorithm  on  sparse  grid,  which  reduced  the  computing  cost  from 
O (nd)  to  ()(n  In4  n). 

(iii) .  High  order  methods  of  SPDES  via  systems  of  forward  backward  stochastic  differ¬ 
ential  equations.  In  this  subproject,  we  overcome  the  problem  of  low  regularity,  by  using  a 
class  of  high  order  numerical  methods.  In  this  approach,  instead  of  solving  SPDEs,  we  solved 
an  equivalent  coupled  system  of  forward  and  backward  stochastic  differential  equations  (BS- 
DEs).  A  special  class  of  Ito- Taylor  formulas  provides  a  high  order  solution  to  the  system  of 
BSDEs.  This  approach,  which  solves  a  system  of  stochastic  ordinary  differential  equations, 
is  fundamentally  different  and  more  powerful  than  existing  methods.  To  the  best  of  our 
knowledge,  our  proposed  approach  represents  a  first  attempt  to  design  high  order  numerical 
algorithms  for  SPDES. 


2  Detailed  description  of  our  accomplishments 

2.1  Quasi  Monte  Carlo  based  Particle  Swarm  Optimization  (PSO) 
method 

A  PSO  algorithm  maintains  a  population  of  M  particles,  and  places  them  in  the  search 
space  of  the  objective  function  which  may  involved  the  solution  of  a  SPDE.  PSO  defines  each 
particle’s  position  as  a  potential  solution  to  the  function  to  be  optimized,  and  then  searches 
for  the  optima  by  updating  its  position  in  every  iterative  step.  Each  particle  is  associated 
with  a  velocity  which  directs  the  flying  of  the  particle  toward  a  new,  presumably  better, 
position/solution.  The  particles  fly  through  the  problem  space  by  following  the  current 
optimum  particles.  In  every  iteration,  each  particle’s  velocity  is  updated  by  following  two 
“best”  values.  The  first  one  is  current  best  solution  it  has  achieved  so  far.  This  value  is 
called  pbest.  Another  “best”  value  that  is  tracked  by  the  particle  swarm  optimizer  is  the 
best  value,  obtained  so  far  by  any  particle  in  the  population.  This  best  value  is  a  global  best 
and  called  gbest. 

After  finding  the  two  best  values,  the  particle  updates  its  velocity  and  positions  with 
following  equations. 

vt+i  =  wvt  +  cpr i  t+i{pbestt  -  at)  +  c2r2jt+1(gbestt  -  at). 
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Convergence  graph  of  f2  from  25-50  iterations 


Convergence  graph  of  f3  from  25-50  iterations 
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Figure  1:  Convergence  graph  for  function  f2 


«t+ i  —  hi  +  vt+1. 

The  key  to  the  success  of  the  PSO  is  choosing  the  initial  particles  effectively.  Traditional 
PSO  methods  use  Monte  Carlo  samples  for  the  initialization  purpose,  which  may  inherit  the 
slow  convergence  feature  of  Monte  Carlo  methods.  In  our  work,  we  use  the  quasi  Monte  Carlo 
methods  to  choose  these  initial  particles  [2],  A  classical  family  of  low-discrepancy  sequences 
are  Halton  sequences  which  are  bases  on  the  radical  inverse  function  defined  as  follows: 


Mn) 


bo  bi 

9 

P  PZ 


+  ...  + 


pm+l  ’ 


(i) 


where  p  is  a  prime  number  and  expansion  of  n  in  base  b  is  given  as  n  =  bo  +  bip  +  ...  +  6mpm, 
with  integers  0  <  bj  <  p. 

Since  Halton  sequence  Xn  in  (0,  l]s  is  defined  as 


Xn  =  {(j)pi (n) ,  </>P2 (n), ...,  (j)Pa (n)), 


(2) 


where  P\,P2,  ■■■,Ps  arc  pairwise  co-primes.  In  practice,  we  always  use  the  first  s  primes  as  the 
bases. 

Comparison  to  other  low-discrepancy  sequences,  Halton  sequences  are  easier  to  imple¬ 
ment.  However,  a  problem  with  Halton  sequence  comes  from  the  correlations  between  the 
radical  inverse  functions  for  different  dimensions.  The  correlations  cause  the  Halton  sequence 
to  have  poor  2-D  projection  for  some  pairing  coordinates.  In  order  to  improve  the  quality 
of  Halton  sequence,  we  use  the  scrambled  Halton  sequence  tobreak  the  cycle  and  correlation 
among  dimensions.  Scrambled  Halton  sequence  can  help  us  to  ignore  the  number  of  points 
and  obtain  good  quality  of  Halton  sequence.  In  Figure  1  we  can  see  that  the  optimal  Halton 
sequence  significantly  speed  up  the  convergence  of  PSO. 
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2.2  Efficient  adaptive  domain  sparse  grid  method  for  SPDES 

In  this  project  we  are  concerned  with  efficient  numerical  solutions  of  the  Zakai  equation,  a 
parabolic  SPDE,  as  follows. 

f  du(t,x )  =  Lu(t,  x)dt  +  Yjk=i  h*k(x)u(t,x)dYt,  ,  . 

\  u(0,x)  =  u0(x),  x  E  Rd, 

where  L  is  the  infinitesimal  generator  associated  with  the  state  process  Xt  given  by 


Lu  = 


1  d2(aa*)iju 

2  dx.  dxj 

i,j  1  0 


E 


dbiU 
dxi  ' 


Zakai  equation  is  a  major  tool  of  solving  nonlinear  filtering  problems.  One  of  the  most 
effective  methods  of  solving  the  Zakai  equation  is  the  split-up  approximation  scheme  where 
the  original  Zakai  equation  is  split  into  a  second  order  deterministic  partial  differential  equa¬ 
tion  in  the  prediction  step,  and  a  degenerate  second  order  stochastic  PDE  in  the  update 
step.  In  the  numerical  simulation  process  of  the  split-up  scheme,  a  prior  PDF  is  obtained 
by  solving  the  deterministic  PDE  at  the  prediction  step;  then  this  prior  PDF  is  updated  at 
the  update  step  following  an  a  posteriori  criterion.  The  strength  of  the  split-up  is  that  it 
is  a  first  order  algorithm  as  opposed  to  most  standard  finite  difference  and  finite  element 
approximations  which  are  half  order,  thus  it  alleviates  the  low  regularity  problem.  However, 
high  dimensionality  and  the  unbounded  domain  are  still  the  main  obstacles  for  the  Zakai 
filter  to  be  an  effective  method  of  solving  nonlinear  filtering  problems. 

In  this  project,  we  constructed  a  fast  split-up  finite  difference  scheme  for  the  Zakai  equa¬ 
tion  using  sparse  grid  and  on  adaptively  chosen  computational  domains.  Specifically,  we  use 
the  importance  sampling  method  to  adaptively  construct  a  bounded  domain  at  each  time 
step  of  the  temporal  discretization.  This  is  a  key  novelty  of  our  algorithm.  Then  we  use 
the  split-up  finite  difference  scheme  to  solve  the  Zakai  equation  on  sparse  grids  over  these 
bounded  adaptive  domains.  The  hierarchical  sparse  grid  method,  which  was  originally  cre¬ 
ated  to  approximate  multi-variable  functions,  uses  only  0(nlogrf-1n)  number  of  grid  points 
instead  of  0(nd)  number  of  grid  points  for  the  standard  full-grid  approximation.  Thus  it  has 
the  potential  of  significantly  alleviating  the  “curse  of  dimensionality”  problem  for  moderately 
high  dimensional  problems. 

To  elaborate  the  basic  idea  of  our  proposed  algorithm,  we  first  create  a  partition  7 Zt  on 

[o,n 

Tit  =  {t„\tn  €  (o.  T],t„  <  tn+i,n  =  0, 1,  ■  ■  •  ,  Nr  -  Mo  =  0,  tNr  =  T} 

and  denote  A tn  =  tn+ 1  —  tn,  n  —  0, 1,  •  •  •  ,  Nt  —  1.  For  n  —  0, 1,  •  •  •  Nt  —  1,  assume  that  un 
is  the  numerical  solution  of  the  Zakai  equation  at  tn.  We  will  use  the  importance  sampling 
method  to  draw  M  sample  points,  denoted  by  ,m,  according  to  the  conditional 

PDF  un  of  the  state  Xn.  Then  we  propagate  each  of  these  samples  from  time  step  tn  to  tn+ 1 
using  the  State  Equation  in  the  nonlinear  filtering  problem  (??)  to  get  M  updated  sample 
points,  denoted  by  {Pn+i}rn=\,  ~  ,m-  Let 

■^n+l  =  [an+\i^n+\\  ^  R 


be  the  smallest  box  domain  containing  all  the  samples  {p™+  \}m=i  an(l  E  =  (E1, . . . ,  Ed)  be 
the  vector  of  marginal  standard  deviations  of  these  samples.  Then,  for  a  user-defined  positive 
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Target  trajectory  and  adaptvie  solution  domains 


our  method  and  the  particle  filter 


constant  A  we  let 

On+l  Qn_|_I  AS 
K+i  =  bn+i  +  AS 

and  choose 

T-^n+l  [®n+lj  ^n+l] 

as  the  solution  domain  for  un+\ ,  which  is  the  approximation  of  the  support  of  un+  \ . 

The  idea  of  the  adaptive  selection  of  the  solution  domain  T>n  is  similar  to  that  in  the 
prediction  step  of  the  particle  filter  method.  As  such  the  region  T>n+i  which  includes  all  the 
updated  samples  is  the  support  where  the  particle  filter  method  builds  the  prior  PDF  of  the 
target  state.  In  our  method,  we  choose  a  confidence  region  surrounding  Vn+ 1  as  our  solution 
domain.  In  this  way,  we  can  use  a  much  smaller  number  of  samples  than  tlie  particle  filter 
method  and  still  maintain  the  accuracy  of  our  algorithm. 

In  the  next  step,  we  use  the  split-up  finite  difference  scheme  on  a  sparse  grid  to  approxi¬ 
mate  the  Zakai  equation  on  each  adapted  domain  T>n. 

We  conducted  several  numerical  experiments  to  demonstrate  our  methodology  in  solving 
nonlinear  filtering  problems.  The  left  subfigure  of  Figure  1  is  the  sketch  of  the  moving  do¬ 
mains  at  several  different  time  steps.  One  can  see  that  the  size  and  location  of  the  computing 
domain  vary  adaptively  at  different  time  steps.  The  right  subfigure  of  Figure  1  are  the  graphs 
of  marginal  conditional  PDFs  at  time  T  =  1  on  the  x  coordinate  using  our  algorithm  for  the 
Zakai  filter  as  well  as  those  obtained  by  the  particle  filter.  From  the  graph  we  can  see  the 
convergence  trend  of  the  particle  filter  toward  the  one  obtained  by  our  algorithm.  When  the 
particle  size  is  160000,  the  solution  generated  by  the  particle  filter  is  very  close  to  our  result 
while  the  results  obtained  by  smaller  sizes  of  particles  are  not  quite  accurate.  In  Table  1  we 
list  the  computing  costs  of  the  particle  filter  with  different  sizes  of  particles,  the  standard 
Zakai  filter,  and  the  Zakai  filter  with  our  algorithm. 

We  can  conclude  hat  our  adaptive  Zakai  filter  algorithm  has  clear  advantages  over  the 
particle  filter:  it  is  more  accurate  with  the  same  computing  cost  and  it  is  far  more  computa¬ 
tionally  efficient  with  the  same  accuracy. 

2.3  Efficient  and  high  order  methods  for  SPDES  through  BSDES 

One  naturally  asks  if  there  is  any  way  to  systematically  construct  higher  order  numerical 
algorithms  for  nonlinear  filtering  problems.  In  this  subproject,  we  answered  this  question  by 


(4) 

(5) 
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Table  1:  CPU  time  comparison 


Method 

CPU  time  (seconds) 

Particle  filter:  20,000  particles 

554 

Particle  filter:  40,000  particles 

1976 

Particle  filter:  80,000  particles 

7703 

Particle  filter:  160,000  particles 

32380 

Standard  Zakai  filter 

270954 

Our  adaptive  Zakai  filter 

587 

taking  a  completely  different  approach  from  existing  methods:  instead  of  solving  the  Zakai 
equation,  we  solve  equivalent  coupled  forward-backward  SDES.  Then  we  construct  a  class 
of  efficient  and  higher  order  numerical  schemes  for  this  coupled  forward-backward  system. 

To  describe  the  basic  idea  of  our  approach,  we  briefly  describe  the  theoretical  results  we 
obtained  in  [?].  These  results  lay  the  foundation  of  our  BSDE  methods  for  solving  SPDES. 
Through  rigorous  stochastic  analysis  we  first  derived  a  system  of  forward  and  backward  SDEs 
(BSDEs)  for  (Xt,  Qt,  Zt) 

I  dXs  =  b(Xa)dt  +  <JsdWs,  Xt  =  x,t<s<T ,  (SDE) 

(dQs  =  ZsdWs-g(Xs)QsdVs,  Qt  =  $(Xt).  (BSDE)  (6) 

Here  Wt  and  Vt  are  two  independent  Brownian  motions.  The  first  equation  in  (6)  is  a 
forward  SDE  while  the  second  equation  is  a  backward  SDE.  One  may  suspect  that  (6)  is 
under- determined  because  it  has  two  equations  but  three  unknown  stochastic  processes  Xt, 
Ut  and  Zt.  However,  under  the  constraint  that  Xt,  Ut  and  Zt  are  adapted  to  some  hltrations 
generated  by  Bt  and  Wt ,  they  can  be  uniquely  determined. 

We  have  a  constructed  a  number  of  numerical  algorithms  to  demonstrate  efficient  of  our 
BSDE  approach  of  solving  SPDES.  Here  we  outline  the  derivation  of  the  Erst  order  scheme 
for  a  general  coupled  system  of  forward-backward  SDEs  [1]: 


'  dXs  =  b(Xs)ds  +  a(Xs)dWs ,  t  <  s  <T, 

dYs  =  +/ (s,  X3,  Ys)ds 

<  +g(s ,  Xs,  Ys)dBs  -  ZsdWs,  t  <  s  <  T,  (7) 

I  Xt  =  x, 

,  UT  =  h(XT). 


For  simplicity  we  assume  that  b  —  0,  a  —  I,  and  g  depends  on  t  and  Y  only.  For  this  simple 
case,  the  two  sided  Ito-Taylor  expansion  for  g  is  given  by 

=  g(o,yo)-f  dy{s,ys)  ■  gd*Ba 

+  fo  g'y(s,  Vs)  ■  zsdws  -  f*  g'y(s,  ys)  ■  fsds 
1  r1  i  U 

~2  Jo  9yy(s,ys- g2)ds+  -  +g'yy(s,ys) 

So  for  a  first  order  algorithm,  the  last  three  terms  are  the  “higher  order”  terms  and  can 
be  dropped  in  the  numerical  algorithm.  Specifically,  if  =  iAt,  i  =  0,  •  •  •  ,  N  is  a  uniform 
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Figure  3:  Comparison  between  the  particle  filter  and  BSDE  filter:  Left:  accuracy  comparison 
in  simulating  the  wells.  Right:  computing  complexity  comparison 


partition  of  [0,  T],  then 

g(tn,Y(tn))  =  g(tn+1,Y(tn+1))  -  g'y(tn+1,Y(tn+1))  ■  Z(tn+1)AWn+1 
+9y(tn+ 1)  Y (tn+i))  ■  g(tn+ 1,  F(fn+i))A5tn+1  +  O(At). 

To  complete  the  construction  of  the  algorithm  we  need  to  eliminate  the  Z  term  in  the 
above  equation.  This  can  be  done  by  taking  the  conditional  expectation  E(-\Jr^s)  where 
T*t  :=  a(x  +  Wr\ t  <  r  <  s)  V  cr(Wt)  V  a{Bt\  0  <  t  <  T).  Then  after  some  manipulations  of 
the  noise  terms  we  derived  the  algorithm  for  the  approximate  solution  Yn  and  Zn  as  follows. 

Yn  =  Yn+ 1  +  AtE[f(tn+i,Yn+i)  —  -g'y(tn+i,Yn+i)  ■  g(tn+i,  Yn+1)\ 

~\~2E[g(tn+ 1?  Yn+ 1)  +  gy(tn+i,  •  g{tnjr\1  Yn+i )]  ABin+1 

AtnZn  =  E[Yn+1AWtn+1]  +  Atn  ■  E[f(tn+i,Yn+i)AWtn+1\ 

+ABtn+ 1  •  (E\g(tn+i,Yn+1)AWtn+1). 

We  proved  in  [1]  that  the  mean  square  error  is  O(At)  for  Yn. 

Theorem  ( Bao ,  Cao,  Meir,  Zhao,  SIAM  J.U.Q.  2015  )  Under  some  regularity  condition 
of  functions  /  and  g,  the  convergence  rate  of  above  algorithm  is  first  order: 

max  E[(Yn  -  Yt  )21  <  C(At)2. 

0<n<N—l  LV  n>  1  K 

In  our  numerical  experiments,  we  apply  our  BSDE  method  to  solve  nonlinear  filtering 
problem.  As  can  be  seen  from  the  left  subfigure  of  Figure  3,  our  BSDE  filter  simulate  the 
switch  much  better  than  the  particle  Liter.  From  the  right  subfigure  of  Figure  3,  we  can  see 
that  with  the  same  root  mean  square  error,  our  BSDE  based  PDF  Liter  requires  far  less  CPU 
time  than  the  particle  Liter.  More  importantly,  our  BSDE  Liter  requires  far  few  steps  than 
that  in  the  particle  Liter  to  complete  the  switch,  which  makes  our  method  far  more  efficient 
than  the  particle  Liter. 
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